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Abstract 



Q^ This article describes the Monte Carlo simulation used to interpret the measurement of the muon-induced neutron flux in the 
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Boulby Underground Laboratory (North Yorkshire, UK), recently performed using a large scintillator veto deployed around the 
ZEPLIN-II WIMP detector Version 8.2 of the GEANT4 toolkit was used after relevant benchmarking and validation of neutron 
production models. In the direct comparison between Monte Carlo and experimental data, we find that the simulation produces 
a 1.8 times higher neutron rate, which we interpret as over-production in lead by GEANT4. The dominance of this material in 
neutron production allows us to estimate the absolute neutron yield in lead as (1.31 + 0.06) x 10"^ neutrons/muon/(g/cm^) for a 
mean muon energy of 260 GeV. Simulated nuclear recoils due to muon-induced neutrons in the ZEPLIN-II target volume (~1 year 
exposure) showed that, although a small rate of events is expected from this source of background in the energy range of interest 
for dark matter searches, no event survives an anti-coincidence cut with the veto. 
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p |1. Introduction 
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Neutrons constitute a very important background for experi- 
ments looking for rare events in deep underground laboratories. 
The vast majority is produced by natural radioactivity in the 
J> rock or the materials used in the detectors: traces of U/Th emit 
a particles which may then interact with light nuclei (Z<30) 
to produce neutrons; another contribution comes from spon- 
taneous fission of heavy elements (mainly ^^^U). The energy 
• range of these neutrons is restricted to a few MeV, and it is 
therefore possible to shield detectors from rock neutrons using 
QQ H-rich materials, which moderate and capture them. Radioac- 
tivity from detector parts can in turn be minimised by an ap- 
propriate choice of building materials and selection of batches. 
Moreover, these local sources of radioactivity can be assessed 
with dedicated measurements (see, for example, Ref. [1] for the 
case of the Boulby Underground Laboratory). Detailed Monte 
Carlo (MC) simulations can then use this information to opti- 
mise the geometry of the passive shielding and to predict this 
contribution to the total background of the experiment, thus re- 
ducing its systematic uncertainties. 

A much smaller contribution comes from muon interactions 
in the rock and surrounding materials, which produce fast neu- 
trons with energies up to the GeV scale. These high energy 
neutrons can easily penetrate through passive shielding (which 
can in fact act as a target for their production) and interact in 
the detectors. In dark matter searches, elastic scattering of high 
energy neutrons produces nuclear recoils within the expected 
energy range of interactions from Weakly Interacting Massive 
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Particles (WIMPs). In double beta decay experiments they con- 
tribute to the y-ray background via inelastic scattering and cap- 
ture; they also mimic neutrino detection in scintillators via in- 
verse beta decay. Active vetoes are thus necessary to help re- 
move this background, by looking at coincidences with high en- 
ergy deposits from the primary muon or the associated cascade. 
In spite of this, these neutrons can travel far from the original 
muon track, interacting in the detector while leaving no telltale 
signature in the anti-coincidence system or the detector itself. 
We show in Section |4] that the rate expected from such events 
in the energy range of interest is quite low for the current gen- 
eration of dark matter experiments based on liquid rare gases, 
which have only a few kilograms of active material. But with 
projects under way to build detectors with hundreds or thou- 
sands of kilograms, the precise knowledge of this neutron flux 
and the ability to model it correctly becomes paramount for the 
design of both the detectors and the associated anti-coincidence 
systems. Moreover, current detectors would also benefit from 
this, as it would allow for a reduction in the systematic errors 
of the expected background. 

Measuring the flux of muon-induced neutrons is not trivial, 
as the requirements are almost those of a low-background ex- 
periment: the detector must be sensitive to neutrons and have a 
large mass, it must be placed in an underground laboratory, and 
it must be stable throughout the duration of the measurement, 
which is typically a few months long. Several experiments have 
tried to measure this flux in the last decades, using either accel- 
erators on the surface \M or dedicated experiments in under- 
ground laboratories fl i, S S H S H El, and several 
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groups have projects for new measurements (see, for example, 
1112.. 13.1 ). But these results are often mutually inconsistent (see 
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lfl4l [Tsi [T^ for a more detailed discussion of available data). 
Moreover, at the time when these experiments were performed 
detailed MC simulations were not possible due to the lack of 
computational power and adequate theoretical models. 

Powerful simulation packages based on advanced theoreti- 
cal models, such as GEANT4 lO andFLUKA Hlllli^l, are 
now available, and although they had some success in explain- 
ing the neutron yield in low-A materials, both fail to explain 
older measurements of neutron production in lead, as was in- 
vestigated in fl?]. Extracting the neutron yield from available 
data is not straightforward, as it requires a detailed description 
of all physical processes involved and of the geometrical setup 
used in the experiment. Moreover, such simulations are not 
trivial either: a good knowledge of the muon flux and energy 
spectrum for the particular underground site is essential, as they 
determine the neutron energy spectrum and yield; the statistical 
efficiency is very low, i.e. a very large number of muons has 
to be simulated in order to produce a statistically significant re- 
sult; and it is all too easy to bias the final result while trying to 
improve this efficiency. 

Recent works have used these simulation packages to try to 
explain published data fl3, or to predict the muon-induced 
neutron background in underground experiments and optimise 
the shielding and veto systems (e.g. i22li23ll ). 

New measurements of total neutron yield in different ma- 
terials, the energy spectrum and lateral distribution regarding 
the primary muon track are therefore urgent, but should be sup- 
ported by detailed MC simulations, preferably using multi-pur- 
pose packages (such as GEANT4 and FLUKA) so that their 
results can be applied to the design of other experiments. 

A new measurement of muon-induced neutrons has recently 
been performed in the Boulby Underground Laboratory (Boulby 
mine. North Yorkshire, UK) |24]. A large mass of lead was 
used as a target that completely surrounded the experimental 
apparatus. The prime goal of this paper is to describe the de- 
tailed Monte Carlo simulations carried out to predict the rate 
of events due to muon-induced neutrons in a scintillator used 
in the aforementioned experiment (Section[3]). Accurate Monte 
Carlo of all physical processes and of the whole setup includ- 
ing trigger conditions and selection cuts, and the comparison 
of the simulation results with the measurements allowed us, for 
the first time, to test theoretical models directly and to make a 
robust prediction of the neutron yield in lead (Section |4|. The 
result for lead is particularly important since i) lead is used as 
a shielding against gamma-rays in high-sensitivity experiments 
for rare event searches, such as direct dark matter, double-beta 
decay and neutrino experiments; ii) previous measurements of 
muon-induced neutron flux in lead reported a higher yield com- 
pared to simulations (carried out, however, by different teams 
not connected to the original experiments). 

We have also performed detailed comparison of our Monte 
Carlo with previous simulations using GEANT4 and FLUKA, 
and investigated the dependence of the neutron yield on the 
models used in simulations (Section |2]i. We show (Section |2]) 
that our simulations agree within a factor of 2 with most avail- 
able data on neutron yields in light targets but cannot explain 
some of the previously reported measurements, in particular in 



heavy targets such as lead. This is probably due to the fact that 
our models were too simplified and did not take into account all 
details of the experimental setup (detector configuration, trig- 
ger, cuts, etc.) that were not publicly available. 

The procedure for muon-induced neutron modelling descri- 
bed in this paper can also be used in similar simulations for any 
other underground experiment for rare event searches. 

2. Model validation 

Processes responsible for neutron production by muons can 
be categorised in four main classes: /) muon capture and ii) 
direct muon spallation result from direct interactions of muons 
with nuclei. In thick targets showers initiated by a muon-indu- 
ced spallation are the dominant source of neutrons: photo- 
and lepto- production, mainly in electromagnetic cascades, ivj 
hadro-production, mainly in hadronic cascades. 

Negative muon capture is only dominant for shallow depths 
(<100 m w.e.), where abundant stopping muons can be cap- 
tured resulting in highly excited isotopes which then emit one 
or more neutrons. In direct muon spallation the muon-nucleus 
interaction may be described through the exchange of a virtual 
photon, resulting in nuclear disintegration. This treatment al- 
lows the use of measured y-N cross sections (real photons) to 
model this process, but breaks down at low energies, when the 
virtuality of the photon becomes comparable to the muon en- 
ergy UM- 

Many particles and mechanisms are involved in neutron pro- 
duction within showers, and so it is not surprising that many 
physics models are required to describe them over different en- 
ergy ranges. In addition, GEANT4 may offer alternative models 
to treat each physical process, which can be selected accord- 
ing to the particular requirements of each simulation (e.g. the 
trade-off between computational efficiency and accuracy is of- 
ten considered). Some of the models available for describing 
interactions of muons and their associated showers were tested 
previously [15] using version 6.2 of the toolkit against other 
simulation results; in this work we use a similar set of models 
with version 8.2: 

• muon-induced spallation: available above 1 GeV muon 
energy; the resulting final states are obtained from pa- 
rameterised hadronic models. 

• photo-production (y inelastic scattering): final states are 
generated by a chiral-invariant phase-space (CHIPS) de- 
cay model below 3 GeV photon energy, while a theoret- 
ical quark-gluon string (QGS) model is used at higher 
energies. 

• hadronic interactions of nucleons, pions and kaons: a 
QGS model is used for high energies and an intra-nuclear 
binary cascade (BiC) for low energies. An older parame- 
terised model (LEP) can be used to cover the intermediate 
range. This is discussed in more detail further ahead. 

• de-excitation ofthe residual nucleus: y and fragment evap- 
oration, fission, Fermi break-up and multi-fragmentation 
(for highly excited nuclei). 
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• neutron interactions are treated using data-driven models 
below 19 MeV. 



produced in inelastic scattering was also considered to prevent 
double counting. 



For electromagnetic interactions the so-called 'Low Energy' 
package was used with production thresholds of a few tens of 
keV for y-rays and ~ 1 MeV for e" and e^ in all materials. Such 
low thresholds are justified by the importance of the photo- 
nuclear process as a significant source of neutrons. 

The importance of the energy scale at which the changeover 
between the BiC and the QGS hadronic models occurs was 
studied. The former aims to describe interactions in the low 
energy range (Ref. |25] recomends its usage below 3 GeV for 
p and n, and 1.5 GeV for pions, pointing out that it should 
work reasonably well up to 10 GeV for p,n). At very low 
energies (below ~7Q MeV, but the real threshold depends on 
atomic mass. A) it reverts to a 'precompound' model, which 
handles the nuclear de-excitation in the pre-equilibrium stage. 
On the other hand, QGS is targeted at high energies (above 
20 GeV) and depends on other models to fragment and de- 
excite the damaged nucleus after the initial interaction. This 
may involve either the precompound model (this association is 
usually known as QGSP) or the CHIPS model (QGSC). There 
is currently no theoretical model available to cover the energy 
region between the BiC and QGS models; the solution com- 
monly adopted is a low-energy parameterisation model (LEP), 
but this is not kinematicaUy correct for a single interaction 12511 . 

We analysed different possibilities to bridge this energy range. 
Firstly the LEP model was avoided altogether for p and n, us- 
ing changeover energies of 6 GeV (first test) and of 10 GeV 
(second test) between the BiC and QGSP models. A third test 
relied on LEP for p and n between 10 and 20 GeV (this is 
very similar to the approach used in the QGSP_BIC_HP ref- 
erence physics list provided by GEANT4 |26]). As for pions 
we never extended the recommended range of the BiC model 
(1 .5 GeV), and thus had to use LEP to bridge to the QGSP min- 
imum energy (note that this is unlike the QGSP_BIC_HP ref- 
erence physics list, where the LEP model is used down to low 
energies). For aU other hadrons, low- and high-energy parame- 
terised models were used. 

These approaches were tested in various materials for differ- 
ent muon energies. We found that the resulting neutron yields 
always agreed within 10%, which is a smaller discrepancy than 
that observed between versions 6.2 and 8.2 of the toolkit. 
GEANT4 offers an alternative model to BiC for the treatment 
of low energy hadronic reactions, the Bertini Cascade |25]: this 
is similarly found not to alter the total yields noticeably [27]. 
We thus opted for a direct changeover at 6 GeV between BiC 
and QGSP for p and n, as it enables a more direct comparison 
with earlier work lll5ll . 



2.1. Total yield in different materials 

We began by benchmarking the neutron yield simulations 
reported in Ref. [15] for different materials using a similar setup. 
This consisted of a beam of mono-energetic /i" incident at the 
centre of a slab of material with thickness 3200 g/cm^; only 
neutrons produced in the central half-length of the slab were 
counted to avoid edge effects. Treatment of secondary neutrons 
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Figure 1 : Variation of tlie neutron yield (per unit muon track length) with the 
initial muon energy for C„H2„ scintillator. Experimental data are taken from 
measurements at various depths between 20 and 5200 m w.e. using the mean 
muon energy for the respective depth. Lines connecting the markers are to 
guide the eye only. 



As most of the neutron yield data available is for organic 
scintillators, we started by studying a Ci()H2o hydrocarbon with 
density p - 0.8 g/cm^. Nevertheless, these results are exten- 
sible to generic hydrocarbons C„H2„, which can also represent 
materials usually used as passive neutron shielding. Figure [T] 
shows the neutron yield as a function of the incident muon en- 
ergy for this material. Statistical uncertainties are comparable 
to the size of the markers. Also shown for comparison are 
the yields obtained with GEANT4 version 6.2 (as reported in 
Ref. [15]) and FLUKA-2008. The neutron yields obtained with 
FLUKA-2008 are in good agreement to what was reported in 
Refs. iQllll using FLUKA-1999, but lower (especially at low 
energies) than the FLUKA-2003 results reported in Ref. lll5ll . It 
is clear that version 8.2 of GEANT4, while still consistent with 
the experimental measurements, produces systematically fewer 
neutrons than version 6.2. At lower energies (< 100 GeV) the 
yield is closer to that predicted by FLUKA, while the agreement 
with version 6.2 of GEANT4 is better at higher energies. 

Figure |2] shows the differential energy spectrum of neutrons 
produced by 280 GeV muons in this material as well as the con- 
tributions from the most important individual processes. The 
photonuclear interaction of y-rays dominates for low energy 
neutrons (<40 MeV), while pion spallation is the biggest con- 
tributor above this energy. The pion absorption process shows 
the expected cut-off right below the pion rest mass subtracted 
by the nucleon binding energy, but unlike what has been re- 
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Figure 2: Differential energy spectrum of neutrons produced in C„H2„ by 
280 GeV muons (in thiclc blaclc). Also shown are the contributions of the most 
important individual processes: photonuclear interaction of y-rays (y — > N), 
neutron inelastic scattering {n — > N), pion spallation {n —> N) and absorption 
(n^abs), and proton (p —> N) and muon (/i — » A') spallation. 



Figure 3: Neutron yield as a function of the average atomic weight for inci- 
dent muons of 280 GeV. Shown for comparison are the results obtained with 
FLUKA-1999 in 1 14], GEANT4 6.2 (and the respective power-law param- 
eterisations) and FLUKA-2008. 



ported in Ref. |23], where the classical GEANT4 nuclear cap- 
ture model was used, does not exhibit a distinct peak before this 
cut-off. In this work we used a new CHIPS-based model, which 
produces a smoother spectrum with more low energy neutrons. 
Muon spallation, responsible for initiating the electromagnetic 
and hadronic showers, clearly contributes very little to direct 
neutron production for muon energies above 100 GeV. Version 
6.2 of GEANT4 yields very similar results for all processes, 
with the notable exception of the photonuclear interaction for 
which it produces more low energy neutrons (<20 MeV). 

We also tested other materials of relevance for low back- 
ground experiments, namely NaCl (which dominates the com- 
position of Boulby rock) and Pb (usually used as passive y-ray 
shielding). For this study we used 280 GeV muons, close to 
the mean muon energy at Boulby (~260 GeV), and kept the 
selection rules mentioned previously. Figure |3] shows the total 
neutron yield as a function of atomic weight (average atomic 
weight for compounds). Results obtained with FLUKA-1999 
H, GEANT4 6.2 [15] and FLUKA-2008 ai-e shown for com- 
parison. For the latter only C„H2„, NaCl and Pb were tested, 
while for the former results for several materials are presented 
along with the respective fits to the power law R = bA^ (R 
is the neutron rate and A is the atomic weight). GEANT4 8.2 
yields are consistent with those obtained with version 6.2, while 
FLUKA-1999 predicts consistently higher rates: ~3Q% higher 
in C„H2„, -50% in NaCl and -80% in Pb. Yields from the new 
FLUKA-2008 are very similar to the ones reported for FLUKA- 
1999 in Ref. |14i] for both C„H2„ and Pb, but much lower for 



NaCl, agreeing with GEANT4 predictions in this case. 

Figure |4] shows the spectra for neutrons produced in each 
of these materials using incident muons of 280 GeV (for NaCl 
and C„H2„) and 260 GeV (for Pb). It is clear that the increase in 
neutron yield seen in materials with large atomic number comes 
mainly from low energy neutrons (<20 MeV). These spectra 
can be compared with the spectra from rock radioactivity (see 
Ref. fP] for the case of Boulby): while neutrons from the rock 
are restricted to a few MeV, muon induced neutrons extend to 
tens and even hundreds of MeV. 

In conclusion, the neutron yields from GEANT4 have not 
changed significantly from version 6.2 to 8.2 for the materials 
and energies of interest for our experimental setup. Differences 
to the FLUKA results are still present and in some cases even 
accentuated, with the latter producing nearly twice as many 
neutrons in lead. As was discussed in Ref. lll5[l this difference 
in the neutron yield may originate from an over-production in 
hadronic cascades by FLUKA: as the emission of fast nuclear 
fragments from highly excited nuclei is not modelled, more en- 
ergy is available for neutron evaporation. 

2.2. Muon-induced spallation 

As mentioned previously, the muon-nucleus interaction can 
be modelled by the exchange of a virtual photon, which allows 
the use of cross-sections parameterised for the case of the real 
photonuclear interaction. For thick targets this is not a dominant 
process for neutron production, and its importance decreases 
with both energy and atomic weight ifTsilTal . Nevertheless we 
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Figure 4: Differential energy spectra for neutrons produced in Pb (thick lines), 
NaCl (dashed lines) and C„H2„ (dotted lines) for 280 GeV incident muons (ex- 
cept for Pb, in which case 260 GeV muons were used). Mean energies of these 
distributions are 8.8 MeV, 23.4 MeV and 65.3 MeV for Pb, NaCl and C„H2„ 
respectively. 



discuss it here in more detail in order to test a new implementa- 
tion based on the CHIPS concept, which has become available 
in recent versions of GEANT4. 

CHIPS (see lEslEgll and references therein) was first avail- 
able in this toolkit as a nuclear de-excitation model, meant to 
be used in conjunction with other models (such as the QGS) 
which handled the initial interaction. Recently, new CHIPS- 
based implementations have become available which include 
the cross-sections of the primary interaction and simulate the 
entire reaction for a number of projectiles including muons. 

As reported in Refs. Uslllil] the older GEANT4 model for 
the simulation of muon-induced spallation (MuNuclear) failed 
to describe the results obtained by the CERN NA55 experi- 
ment which claimed to have measured fast neutron pro- 
duction in graphite, copper and lead using a beam of 190 GeV 
muons. Neutron detectors with a reported threshold of ~10 
MeV were placed at 45°, 90° and 135° regarding the initial 
muon beam. However, the thin target assumption in this exper- 
iment has been questioned We investigated whether the 
newly available models could explain the discrepancies found 
previously. Four simulations were conducted: a) using only the 
previous GEANT4 muon spallation process; b) using the com- 
plete set of physics processes described previously; c) using 
only the new CHIPS model for the simulation of muon pho- 
tonuclear interactions; and d) using the complete set of physics 
processes, but selecting CHIPS versions whenever available (na- 
mely for all lepto-nuclear interactions, capture of negatively 
charged particles at rest, and nuclear de-excitation). Note that 



a) and c) should produce results similar to b) and d), respec- 
tively, if the thin target assumption is correct. 

Figure|5]shows the differential cross-sections obtained from 
considering only the primary muon interaction (left) or all the 
physics processes (right). Results with the standard physics 
processes (a) and b)) reproduce those obtained with older ver- 
sions II15L 12111 . and fail to explain the experimental results. The 
new CHIPS model for muon spallation produces a dramatic in- 
crease for graphite, showing a reasonable agreement with the 
data; using the complete set of physics processes does not chan- 
ge this result. As for copper, the primary process by itself al- 
ready produces more neutrons than the older one, but in this 
case it is only the use of the complete physics list that brings 
the simulation result closer to the measured yields, showing that 
this could not be considered as a thin target. 

The very large discrepancy observed for lead remains unex- 
plained, with the CHIPS-based models yielding approximately 
the same result as the older versions, about one order of magni- 
tude lower than the measurement. Considering a lower energy 
threshold for the neutron detectors helps to reduce this differ- 
ence, but a value as low as 3 MeV is required for a reasonable 
agreement. In any case the shape of the distribution still dis- 
agrees with the measurements. 

Figure|6]compares the neutron energy spectra using physics 
lists with complete sets of processes (cases b) and d)) for each 
of these materials. Results from both simulations are more sim- 
ilar in the case of lead, but as already indicated from Figure |5] 
completely fail in reproducing the experimental results for all 
three angles. CHIPS produces a better agreement for graphite, 
but is still lower for 45° at high energies. The differences be- 
tween simulation and data are not as evident in copper as in 
lead, but both models still clearly fail in describing their be- 
haviour. 

From the analysis of Figures |5] and |6] it is clear that none 
of the tested models is able to describe the whole set of results 
from the NA55 experiment. However, this disagreement should 
be taken with caution: it is difficult to do such simulations retro- 
spectively based on the little information available in the litera- 
ture. This reafirms the need for new experiments accompanied 
by dedicated and detailed MC simulations of the experimental 
setup and selection criteria. 

2.3. Muon transport 

We also tested muon transport in GEANT4 against FLUKA- 
2007 and the muon propagation code MUSIC ll30ll . For the 
case of GEANT4 the two models described in l2.2l to handle di- 
rect muon spallation (MuNuclear and CHIPS) were used, but 
no other change was done to the physics list. 

First, one million muons with 100 GeV and 1000 GeV were 
propagated through slabs of different materials with a thickness 
equivalent to 10 m w.e. Figure |7] shows the energy distribu- 
tions of 1000 GeV muons upon exiting a slab of water Both 
GEANT4 models agree well with MUSIC and FLUKA, with 
the only exception of a factor of 2 enhancement in the num- 
ber of muons below 50 GeV (5% of the initial energy) for the 
CHIPS model. We would hke to point out, however, that, since 
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Figure 5: Differential cross-section of neutron production in (curves from bottom to top) graphite, copper and lead by 190 GeV muons (for a 10 MeV threshold). 
Plot on the left shows GEANT4 simulations considering only the direct muon spallation interaction - cases a) and c) described in the text - while the one on the 
right includes the complete physics lists - cases b) and d). Solid lines show standard process results; dotted lines the new CHIPS-based processes; markers show 
the NA55 results. 



electromagnetic processes, such as ionisation, bremsstrahlung 
and pair production dominate in muon energy losses, the effect 
of muon inelastic scattering (or hadroproduction) can hardly be 
seen with this approach. 

We have also checked the reliability of GEANT4 to prop- 
agate muons through a large thickness of matter This is rel- 
evant for the calculation of muon fluxes and spectra at large 
depths underground, and to predict muon event rates from high- 
energy neutrino interactions. One million muons with 2 TeV 
energy have been transported through 3 km w.e. of standard 
rock (Z = 11, A = 22, density = 2.65 g/cm^). Figure[8]shows the 
resulting energy distributions from the three codes: GEANT4, 
MUSIC and FLUKA. The two GEANT4 models predict very 
similar energy spectra and only the results using MuNuclear 
are shown in the figure. Here again the differences between 
the two models for muon inelastic scattering are hidden in the 
much larger energy losses due to electromagnetic interactions. 
The probability for a 2 TeV muon to survive after 3 km w.e. 
in standard rock is: 0.730 (MUSIC), 0.750 (GEANT4) and 
0.741 (FLUKA), with a typical statistical error of 0.1%. Mean 
muon energy after propagation is: 261 GeV (MUSIC), 256 GeV 
(GEANT4), 273 GeV (FLUKA), with a statistical uncertainty 
of less than 1 GeV. 

The difference between the MuNuclear and CHIPS models 
can be seen in Figure|9] where the energy distributions are plot- 
ted after muon propagation taking into account only the muon 
inelastic scattering process. Muons with 2 TeV energy have 
been transported through 3 km of water and their energy spectra 
have been recorded. CHIPS model predicts larger muon energy 



loss, resulting in the larger number of muons with smaller ener- 
gies compared to the MuNuclear model or MUSIC. A large en- 
hancement in the number of muons is seen in CHIPS below 100 
GeV (5% of the initial muon energy), an indicator that this pro- 
cess is highly suppressed for low-energy muons in this model. 

3. Simulation of the experiment 

3.1. Experimental setup 

The measurement of the muon-induced neutron flux in the 
Boulby Underground Laboratory (1070 m deep, 2850 w.e. [2?]) 
was performed using a liquid scintillator detector which served 
as a veto system for the ZEPLIN-II dark matter detector ll3lll . It 
consists of a large metal container surrounding the lower half of 
the ZEPLIN-II detector, filled with 0.73 tonnes of liquid scin- 
tillator with a known composition and density. It was used to 
detect both primary high-energy muons (used as trigger) and 
secondary (low-energy) y-rays from neutron capture using de- 
layed coincidences. Independent electronics and data acquisi- 
tion systems were installed for this measurement, running in 
parallel with the ZEPLIN-II experiment (see Figure [TOli. 

ZEPLIN-II and its anti-coincidence system are surrounded 
by a shielding structure made of lead which was designed to 
attenuate y-rays emitted from the rock walls. This lead 'castle' 
has a thickness of 15 cm on the top section and 22.5 cm on the 
side walls and floor. Weighing in at approximately 50 tonnes, 
it creates an excellent target for the production of neutrons by 
muons. Between the top of the castle and the ZEPLIN-II de- 
tector, Gd-impregnated wax (0.2% Gd by weight) and a thick 
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Figure 6: Energy spectra of neutrons produced in (top to bottom) graphite, copper and lead for 135°, 90° and 45° (1st, 2nd and 3rd columns, respectively). Thick 
lines show simulation results usingCHIPS models - case d) described in the text; thinner lines using the models described in section[2]- case b). Data points were 
taken from the NA55 publication 
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Figure 7: Energy spectrum of 1000 GeV muons after crossing a 10 m thick slab 
of water. 



polypropylene sheet were used for the shielding of rock neu- 
trons. Several layers of polypropylene (interleaved with thin 
layers of Gd-loaded wax) were also used inside the lead castle, 
above the veto system. Furthermore, the inner surface of the 
veto vessel was painted with a mixture containing Gd salt. The 
purpose of using Gd around the detector was to increase the 
number of neutron captures, enhancing the probability of neu- 
trons produced by radioactivity in detector components being 
detected by the veto or the ZEPLIN-II detector itself. 

Further details on this experiment, including data acquisi- 
tion and analysis procedures, can be found in Ref. |24] and 
references therein. 

3.2. Simulation details 

The simulation of this experiment can be divided in two 
stages. In the first, the MUSUN lll411 code was used to sample 
muons according to their energy spectrum and angular distri- 
bution at the Boulby Underground Laboratory (obtained with 
the code MUSIC 13011 propagating a spectrum of atmospheric 
muons from the surface through Boulby rock). The muon flux 
thus obtained is normalised to the one measured in this exper- 
iment, which agrees within 10% with previous measurements 
1 3211 . This results in a mean energy of ~260 GeV for Boulby 
underground muons. These primaries were sampled on the sur- 
face of a parallelepiped encompassing the experimental hall, 
with dimensions chosen to give a clearance of 10 m, 7 m and 
4 m to the laboratory walls (top, sides and bottom, respectively). 
Two million of these muons were generated and their proper- 
ties (energy, momentum, position and charge sign) stored for 
the second stage. 



Figure 8: Energy spectrum of 2 TeV muons after crossing 3 km w.e. of stan- 
dard rock (Z=ll, A=22, p=2.65 g/cm^). Mean muon energies obtained from 
MUSIC, GEANT4 and FLUKA ai'e 261, 256 and 273 GeV, respectively, with 
survival probabilities of 0.730, 0.750 and 0.741. 



A detailed GEANT4 simulation of the experiment was de- 
veloped for the second stage using version 8.2 of this toolkit. 
The experimental hall, lead castle, veto detector and the en- 
tire ZEPLIN-II detector were included, as well as the neutron 
shield (with Gd-loaded materials). A thin layer of Gd was also 
applied to the inner surfaces of the veto, in order to mimic the 
Gd-loaded paint. FigurefTOl shows a cross-sectional view of the 
final geometry obtained with GEANT4. 

Muons sampled in the first stage are the input to this simula- 
tion. The primary muons and all secondary particles created in 
the muon-induced cascades are propagated inside this geome- 
try using the set of physics models described in Section|2l Each 
energy deposit inside the veto is recorded in a time histogram, 
which can be directly compared with the experimental results. 
Energy deposits inside the ZEPLIN-II active volume are also 
recorded to allow coincidence studies. For neutron captures, 
the capture material and isotope are stored, together with infor- 
mation about the production of the captured neutron: material, 
parent particle and production process. For the case of inelas- 
tically scattered neutrons we postulate that the highest energy 
neutron in the final state keeps the identity of the incident neu- 
tron. Capture times are also stored, so that they can be cross- 
referenced to energy deposits inside the veto. 

About 40% of the experimental data was recorded while the 
top section of the lead castle (and the Gd-loaded wax, C in Fig- 
urefTOli and the horizontal polypropylene sheet on top of the de- 
tector were not in place - these are called 'roof-off ' runs, as op- 
posed to the usual 'roof-on' runs performed with the complete 
shielding. This modification to the geometry was also included 
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Figure 10: Vertical cut of the geometry model used in the GEANT4 simulation: A - ZEPLIN-II detector, B - liquid scintillator detector (veto), C - Gd-loaded 
wax, D - lead castle, E - polypropylene sheets which make up the passive neutron shielding (vertical slabs are interleaved with Gd-loaded resin). Details of the 
ZEPLIN-II detector were removed from this figure for simplicity. 



in the simulation for a similar fraction of the processed events. 
In total some 120 million muons were simulated, correspond- 
ing to a total live time of 960 days, about 4.7 times longer than 
the experimental exposure, with each muon generated in the 
first stage being transported through the rock and setup about 
60 times using different sets of random numbers. 

5.5. Validation 

The detector was calibrated in energy with a ^'^Co source 
(at the beginning and end of the run); its sensitivity to neu- 
tron captures was tested with a 0.1 GBq (a activity) Am-Be 
source at the beginning of the run. Both tests were simulated 
with GEANT4 using the same geometry and set of physics 
models described above, with the only difference being the pri- 
mary particles passed to the simulation: for the energy calibra- 
tion y-rays were emitted with energies from either of the two 
«"Co lines (1.173 MeV and 1.333 MeV); for the Am-Be run 
the SOURCES4A fs?] code was used to generate the energy 
spectrum of the initial neutrons. 

Comparison of these simulations with experimental data 
shows a very good agreement (particularly in terms of time de- 
lay distribution of y-rays from neutron captures) |24], thus val- 
idating the geometry implementation and physics models used. 



Note, however, that the neutron calibration is not ideal for eval- 
uating the detector response to muon-induced neutrons: Am- 
Be neutrons have, on average, lower energies and are created 
at a fixed location (just below the castle roof), whereas neu- 
trons originated by muons have higher energies and can be cre- 
ated anywhere. There are also other differences, due to exper- 
imental limitations: the triggers used are necessarily different 
(neutron-induced proton recoils instead of high-energy deposits 
from muons) and the fact that the acquisition system was run- 
ning with almost 100% dead-time during calibration. All these 
problems are obviously not particular to this experiment, but 
rather apply to any similar experiment that uses neutron sources 
for their calibration. 

Nevertheless, this calibration allows us to show that this ex- 
periment is sensitive to neutrons and that our simulation is ac- 
curate enough to simulate neutron transport and detection in the 
MeV and sub-MeV energy ranges. 

4. Results 

Some results from this simulation were already analysed in 
Ref. [24], where it was used as an aid in the interpretation of 
experimental results. 
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Figure 9: Energy spectrum of 2 TeV muons after crossing 3 km of water. Only 
tlie muon inelastic scattering process was considered. The GEANT4 CHIPS 
model shows a large increase in the number of muons below 100 GeV, an indi- 
cation that the model may be unsuited for this energy scale. 



Comparison of the simulated muon rate with the experi- 
mental value of 52.9 + 0.5 per day allows the estimation of the 
muon flux at Boulby as (3.79 + 0.04 (stat) + 0.1 1 (syst)) x 10"^ 
cm^^s with the systematic error coming from an uncertainty 
in the energy scale, discussed in detail in Ref. (l^. The verti- 
cal depth obtained from this measurement is 2850 + 20 m w.e., 
with the assumption of a flat relief on the surface. Detailed de- 
scriptions of the procedures to estimate the muon flux and depth 
can be found in Ref. js^l- 

As the light collection in the veto was not simulated, the 
energy deposition was smeared using a Gaussian distribution 
with a standard deviation, cr, given by the equation cr/E = 
^fcT+'pjE [34]. The parameters a and ji were determined by 
comparing the smeared spectrum with the measurements, and 
the result was used for the energy calibration of the experimen- 
tal data and the definition of the threshold for muons (20 MeV) 
and delayed y-rays (0.55 MeV). 

The shape of the time delay distributions (of y-rays from 
neutron captures) for 'roof-on' and 'roof-off' simulations was 
compared with experimental data (using the combined data from 
both experimental runs), showing a good agreement in the in- 
terval 40 - 190 //s. The first 20 /is were not considered due to 
the high level of after-pulsing in the experimental data follow- 
ing the muon signal, while the period from 20 to 40 //s was re- 
moved because of the high probabihty of detecting y-rays from 
neutron captures in gadolinium, which is difficult to simulate 
accurately due to uncertainties in its distribution in the wax, 
resin and paint materials used. 

The agreement in the shape of the simulated and measured 



time delay distributions from 40 - 190 /is also enables a direct 
comparison of the neutron rates (defined as the average num- 
ber of detected neutrons per detected muon), by counting the 
y-rays from captured neutrons inside this interval in both distri- 
butions. The average rate obtained from the experimental data 
is 0.079 + 0.003 neutrons/muon (0.084 + 0.004 neuti'ons/muon 
for 'roof-on', 0.072 ± 0.005 neutrons/muon for 'roof-off'). The 
average rate from the simulation is 0.143 + 0.002 + 0.009 neu- 
trons/muon (0.143 + 0.002 neutrons/muon for 'roof-on', 
0.143 + 0.003 neutrons/muon for 'roof-off'). The systematic 
uncertainty comes mainly from the already mentioned uncer- 
tainty in the energy scale. There is a factor of 1.8 difference be- 
tween simulated and experimental values, with GEANT4 pre- 
dicting a higher neutron rate than measured. Given the agree- 
ment in the shape of the time delay distribution in both the Am- 
Be neutron calibration and the background runs, it is unlikely 
that this difference results from a problem in the geometry setup 
of the simulation or the capture and low energy transport mod- 
els of neutrons (these models were also tested in Ref. Issll and 
found to be in good agreement with MCNPX ll36ll ). Therefore, 
the most probable explanation for the excess of muon-induced 
neutrons lies in the GEANT4 models involved in neutron pro- 
duction. 

Figure [TT] shows the relative contribution of each neutron 
multiplicity to the overall number of detected y-rays for both 
experimental and simulated data. It is clear that GEANT4 pro- 
duces an excess of neutrons for all multiplicities up to m = 
12. Constraints in the data acquisition system effectively limit 
the number of recorded pulses to 16, but analysis of simulated 
events with larger multiplicities show that the fraction of missed 
neutrons is at the level of a few percent, and thus does not con- 
tribute significantly to the observed difference in the rates. 

Table [T] shows the contributions of the most relevant ma- 
terials for neutron production. In the 'detected' columns only 
neutrons with y-ray hits above threshold and in the time inter- 
val 40 - 190 fis were used, while 'all' includes all neutrons that 
ended up being captured in the simulation (although very small, 
there is a finite probability that a neutron leaves the geometrical 
limits of the simulation without being captured). As expected, 
nearly all neutrons created in this experiment come from rock, 
but only fewer than 10% are detected, attesting the effectiveness 
of the shield. Lead is by far the biggest producer of detected 
neutrons, being responsible for ~90%. Even in 'roof-off' runs, 
when there is a 12% reduction in Pb mass, its contribution to the 
number of detected neutrons decreases by less than 1%. This 
can be due to two reasons. First, the presence of a large volume 
of Gd-loaded wax in the roof section enhances the moderation 
and capture of neutrons in this region, further away from the 
veto than the H-rich materials. Moreover, the chosen time win- 
dow deliberately excludes the period in which the captures in 
Gd dominate. Scintillator itself does not give a noticeable con- 
tribution to neutron production. 

Lead contributes 90% to the neutron production, making 
uncertainties in other materials and elements of the geometry 
relatively unimportant. The increased neutron rate in the sim- 
ulation therefore implies that GEANT4 overestimates the neu- 
tron production in this material. We note that FLUKA predicts a 
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Figure 1 1 : Relative contribution of each multiplicity to the total number of de- 
tected y-rays from neutron captures in experimental (solid black) and simulated 
data (shaded). 

Table 1 : Contributions of different materials to neutron production for captured 
neutrons. 'Detected' means neutrons which deposited at least 0.55 MeV in the 
veto (after gaussian smearing) during the interval 40 - 190 /is in events with a 
trigger signal (muon energy deposition) of at least 20 MeV; 'all' refers to all 
captured neutrons. 'Roof-on' and 'roof-off' refer to, respectively, runs with 
the complete shielding and runs without the roof sections of the Pb castle and 
neutron passive shielding. 



Material 


detected 


all 




roof-on roof-off 




Lead 


91.4% 90.6% 


0.6% 


Rock 


6.1% 8.0% 


99.4% 


Copper 


2.5% 1.4% 





higher rate, again 80% in excess of GEANT4. This also contra- 
dicts the NA55 results, as discussed in Section l272l Moreover, 
the alternative (CHlPS-based) models yield a ~ 1 .6x higher rate, 
and thus go in the wrong direction regarding these experimental 
results. 

From the simulation we find that only 42% of neutrons are 
detected in the time window 40 - 190 yus, and can thus esti- 
mate the experimental rate for an infinite time window to be 
0.188 + 0.005 neutrons/muon. Considering that most neutrons 
are created in lead and that we have confidence in the simulated 
geometry and the GEANT4 models that treat neutron modera- 
tion and capture, and assuming that the ratio of simulated-to- 
measured neutron rates is the same as for raw neutron yields 
in the material, it is possible to estimate the absolute neutron 
yield of 260 GeV muons in lead. Using the neutron yield ob- 
tained in Section im for 260 GeV neutrons (2.37 x 10"^ neu- 
trons/muon/(g/cm^)) and the 1.8 ratio between the simulated 



and measured rates, we estimate an actual neutron yield in lead 
of (1.31 + 0.06) X lO"-' neuti-ons/muon/(g/cm^), which is 3.2 
times smaller than expected from FLUKA-2008. 

In Table|2]the contribution of different elements to the cap- 
ture of neutrons is shown for the two situations described above 
('detected' and 'all'). As expected, hydrogen is responsible for 
almost two thirds of the detected neutrons, again with little dif- 
ference between the 'roof-on' and 'roof-off' situations. Scin- 
tillator is obviously the main contributor to detected hydrogen 
captures, with 85% of the neutrons being captured in this ma- 
terial, while polypropylene is responsible for only 15%. The 
small contribution from Gd is again due to the choice of the 
time interval for detected pulses. Na and CI (the main con- 
stituents of the Boulby rock) are responsible for almost all the 
neutron captures (99.7%), but have a small contribution to the 
detected y-rays even when the roof of the lead castle was not in 
place, an effect of the reduced solid angle and the selection cuts 
used in the experimental data (threshold and time interval). 

Table 2: Contributions of different elements to neutron captures. 'Detected' 
means neutrons which deposited at least 0.55 MeV in the veto (after gaussian 
smearing) during the interval 40 //s to 190 /iS in events with a trigger signal 
(muon energy deposition) of at least 20 MeV; 'all' refers to all captured neu- 
trons. 'Roof-on' and 'roof-off' refer to, respectively, runs with the complete 
shielding and runs without the roof sections of the lead castle and neutron pas- 
sive shielding. 



Element 


detected 
roof-on roof-off 


all 


H 


65.7% 


64.4% 


0.1% 


Gd 


12.5% 


10.0% 


0.1% 


CI 


12.1% 


14.0% 


94.5% 


Cu 


4.0% 


4.4% 




Fe 


4.0% 


3.4% 




Na 


0.3% 


0.5% 


5.2% 


others 


1.4% 


3.3% 


0.1% 



The contribution of muon-induced neutrons to the back- 
ground of the ZEPLlN-11 detector was also estimated using 
this simulation. For this study 55.5 million muons were simu- 
lated (corresponding to 1 .23 years exposure) using the 'roof-on' 
geometry, and energy depositions in the xenon were recorded 
along with the hits in the veto to allow anti-coincidence stud- 
ies. Moreover, nuclear recoils (NR) were recorded separately 
and normalised by a conversion factor to obtain the equivalent 
energy for electromagnetic interactions (Epf): E^r = Eeg/Q.36 

The results of this study are summarised in Table |3] where 
the contribution of different types of NR events to the back- 
ground is shown, using thresholds and energy ranges of interest 
for dark matter experiments (namely the interval 2-20 keV 
where most WIMP interactions are expected to occur). The ta- 
ble lists both 'pure' and 'mixed' recoil events, the latter repre- 
senting events that involve both NR and electromagnetic energy 
depositions. 

The predicted number of events involving at least one nu- 
clear recoil (and with any energy deposition) in the target vol- 
ume of ZEPLlN-11 is 57.7 + 6.9 per year, of which 13 + 3.3 
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Table 3: Muon-induced neutron background (NR events per year) in the 
ZEPLIN-II (31 kg) target volume for different detection thresholds and energy 
intervals. 



Event type 



Eee (keV) Hits/year 



Mixed NR events 


> 





57.7 ± 0.9 


Pure NR events 


> 





13 + 3.3 






9 

z. 


J.J ± l.O 




2- 


20 


3.3 + 1.6 


Pure single NR events 


> 





5.7 + 2.2 




2- 


20 


0.8+0.8 


Anti-coincidence with veto 


> 





< 1.9 (90%C.L.) 



have no coincident electromagnetic deposition. 3.3 ± 1.6 of 
these events are expected to fall in the interesting energy range 
for WIMP search, but this number can be further reduced to 
0.8 ± 0.8 if we consider that the detector has enough spatial 
resolution to resolve the locations of multiple neutron elastic 
interactions and remove these events. Nevertheless, when we 
introduce an anti-coincidence cut with the veto signal (even us- 
ing a threshold as high as 1 MeV) no NR event survives. If 
only the inner volume of ZEPLlN-11 is used in the data analy- 
sis (7.2 kg) these neutron rates are reduced by approximately a 
factor of 4. 



Data from the ZEPLlN-11 science run Bill , with a total ex- 
posure of 225 kgxdays, was re-analysed to search for neutron 
recoils in coincidence with a large muon signal in the veto, but 
none was found. 

These results show that the muon-induced background is 
not a severe threat for the current generation of detectors, with 
target masses up to a few tens of kilograms, but may become a 
problem for larger scale targets. 

5. Conclusions 

Version 8.2 of the GEANT4 toolkit was used to create a de- 
tailed Monte Carlo simulation of the first experiment to measure 
the muon-induced neutron background at the Boulby Under- 
ground Laboratory, which used the 0.73 tonne liquid scintilla- 
tor veto from the ZEPLlN-1 and ZEPLlN-11 dark matter search 
experiments. 

First, the neutron yields of mono-energetic muons in dif- 
ferent materials of interest to this experiment were compared 
with results obtained with older versions of this toolkit and 
other simulation packages, to benchmark recent modifications 
to the relevant models. Yields are consistently lower than the 
FLUKA predictions for all three materials tested and 280 GeV 
muons, the differences being 30% for C„H2„, 50% for NaCl 
and 80% for Pb, but are still in good agreement with version 
6.2 of GEANT4. Newly available CHlPS-based models were 
also tested while attempting to describe results from the NA55 
experiment, which tried to measure the neutron production in 
thin targets of different materials. These new models clearly 



produce higher yields for all materials and detection angles, but 
are still far from explaining the controversial NA55 results. 

A direct comparison of the simulated and measured neutron 
capture rates in the 40-190 yus window after the trigger reveals 
that the simulation predicts 1.8 times more neutrons than mea- 
sured (well beyond statistical and systematical errors), with a 
similar contribution from mutiplicities up to 12. Given that Pb 
is responsible for ~9Q% of the detected neutrons, this implies 
that the GEANT4 models over-predict the neutron production 
in this material, contrary to what was expected from the com- 
parison with FLUKA predictions and the NA55 results. The 
neutron yield in lead was estimated to be (1.31 + 0.06) x 10"^ 
neutrons/muon/ (g/cm^) . 

Finally, the contribution of muon-induced neutrons to the 
ZEPLlN-11 background was estimated. The expected single 
nuclear recoil rate is less than one event per year, but an anti- 
coincidence with the veto effectively removes all nuclear recoil 
events, confirming that this is not a meaningful source of back- 
ground for a detector of this size. 
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